function [phiest, phiml, sghest, sghml, clayest] = HydrateInversion(data,traindata,ncomp)

% Input: data = Inverted Vp Vs Rho
% Input: ncomp = number of component in the mixture (try 2 or 4)
% Output: phi, sgh = estimated phi and sgh


%% Data
ns=size(data(:,1),1); % number of samples

%% Joint Gaussian Mixture
% joint random variable (5D)
mdjoint=traindata;
nmd=size(mdjoint,2); % number of model + data variables
nd=size(data,2);
nm=nmd-nd;

% components
f_post = kmeans(traindata(:,nm+1:end),ncomp);
% f_post = ones(size(traindata,1),1);
% f_post(mdjoint(:,1)<=mean(mdjoint(:,1)) & mdjoint(:,2)>=mean(mdjoint(:,2)))=2;
% f_post(mdjoint(:,1)>mean(mdjoint(:,1)) & mdjoint(:,2)<mean(mdjoint(:,2)))=3;
% f_post(mdjoint(:,1)>mean(mdjoint(:,1)) & mdjoint(:,2)>=mean(mdjoint(:,2)))=4;

% joint distribution 
mu_joint=zeros(ncomp,nmd);
sigma_joint=zeros(nmd,nmd,ncomp);
for i=1:ncomp
    mu_joint(i,:)=mean(mdjoint(f_post==i,:));
    sigma_joint(:,:,i)=cov(mdjoint(f_post==i,:));
end
p=1/ncomp*ones(ncomp,1);

%% Conditional Gaussian mixtures
mu_cond=zeros(ns,ncomp,nm);
sigma_cond=zeros(nm,nm,ncomp);
pf_post=zeros(ns,ncomp);
f=zeros(ns,1);
phiml=zeros(ns,1);
phiest=zeros(ns,1);
sghml=zeros(ns,1);
sghest=zeros(ns,1);
clayest=zeros(ns,1);

for i=1:ncomp
    sigma_cond(:,:,i)=sigma_joint(1:nm,1:nm,i)-sigma_joint(1:nm,nm+1:end,i)*(sigma_joint(nm+1:end,nm+1:end,i)\sigma_joint(nm+1:end,1:nm,i));
end

for j=1:ns
    for i=1:ncomp
        mu_cond(j,i,:)=mu_joint(i,1:nm)'+sigma_joint(1:nm,nm+1:end,i)*(sigma_joint(nm+1:end,nm+1:end,i)\(data(j,:)-mu_joint(i,nm+1:end))');
        lambda_sum=0;
        for k=1:ncomp
            lambda_sum=lambda_sum+p(k)*mvnpdf(data(j,:),mu_joint(k,nm+1:end),sigma_joint(nm+1:end,nm+1:end,k));
        end
        pf_post(j,i)=p(i)*mvnpdf(data(j,:),mu_joint(i,nm+1:end),sigma_joint(nm+1:end,nm+1:end,i))./lambda_sum;
    end
    [~, f(j)]=max(pf_post(j,:));
    phiml(j)=mu_cond(j,f(j),1);
    sghml(j)=mu_cond(j,f(j),2);
    for i=1:ncomp
        phiest(j)=phiest(j)+pf_post(j,i).*mu_cond(j,i,1);
        sghest(j)=sghest(j)+pf_post(j,i).*mu_cond(j,i,2);
%         clayest(j)=clayest(j)+pf_post(j,i).*mu_cond(j,i,3);
    end
end